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Abstract 

We present EKHARA Monte Carlo event generator of reactions e + e~ — > 
e + e~7r° and e + e~ — > e + e~7r + 7r~ . The newly added channel (e + e~ — > e + e~7r°) 
is important for 7*7* physics and can be used for the pion transition form 
factor studies at meson factories. 
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and Empirical Models and Theories 
Nature of problem: 

The first version of EKHARA [1,2] was developed to simulate background for 
the pion form factor measurement at meson factories coming from the process 
e + e~ — > e + e~ir + ir~ . The newly added channel e + e~ — > e + e~7r° will help in the 
pion transition from factor studies at meson factories [3]. 
Solution method: 

Events consisting of the momenta of the outgoing particles are generated by 
Monte Carlo methods. The generated events are distributed accordingly to the 
theoretical cross section. For the e + e~ — > e + e~7r° mode the Monte Carlo sampling 
developed in [4] was adopted. 
Restrictions: 

In order to compile the code, the FORTRAN 77 compiler should support 
quadruple precision numbers. 
Unusual features: 

Calculations are carried in quadruple precision, in order to avoid numerical 
cancellations in e + e~ — > e + e~ir + Tr~ mode. 

Running time: 

Depends on the requested mode and applied kinematic cuts. Example: on Intel 
Core2 Quad CPU Q6600 @ 2.40GHz, using only one thread, 

• 10 5 unweighted e + e~7r° events are generated in 78 seconds (no cuts), 

• 10 5 weighted e + e~7r + 7r~ events are generated in 38 seconds (with cuts [2]). 
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LONG WRITE-UP 



1. Introduction 

1.1. The physics case 

The first version of EKHARA MC generator 0, Q (EKHARA ver. 1) 
was developed to simulate the reaction e + e~ — > 7r + 7r~e + e~. This process is a 
background [3J to the e + e~ — > 7r + 7r~7(7) cross section measurement at meson 
factories when only charged pions are observed [4], |5j. Its proper simulation 
is relevant for the precise extraction of the cross section a(e + e~ — > 7c + tt~) 
and the charged pion form factor using the radiative return method. The de- 
scription of a physics content (matrix elements, modelling of the pion-photon 
interactions, etc.) in EKHARA ver. 1 was given in 0], however, the compu- 
tational issues remained unpublished. In view of the accepted prolongation 
of the experiment, namely, the KLOE-2 project ja , and the planned radia- 
tive return program of the BES-III experiment 0,lS], this gap has to be filled 
and the description of the e + e~ — > 7r + 7r~e + e~ mode of EKHARA generator 
is one of the aims of this paper. EKHARA ver. 1 is also a convenient base 
for an inclusion of other channels of the inelastic e + e~ scattering. 

The KLOE-2 project assumes an installation of special tagging devices (sf, 
which will allow [6] to study, among other, the e + e~ —> n°e + e~ process, 
aiming at the following measurements: 



the two-photon decay width of 7r 



the transition form factor F 7r o 7 * 7 * (m^. , qf, q%) at space-like photon mo- 
mentum transfers qf, q\. 



These measurements are of a significant importance [6] and the knowledge 
of F 7r o 7 * 7 * for high photon virtualities is supposed to help in the reduction 
of an error in the calculation of the light-by-light contributions to the muon 
anomalous magnetic moment. A reliable Monte Carlo generator for this 
channel, based on the modern knowledge of the tt° — 7* — 7* transition form 
factor is an indispensable tool for such studies. A description of this new 
mode, which is distributed within EKHARA ver. 2, is the second aim of this 
paper. 
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1.2. Basic functionalities of EKHARA ver. 2 

• e + e~ — > e + e~7T + 7T~ 

- generates weighted events; 

- fills the histograms; 

- gives the integrated cross section within cuts. 

• e + e~ — > e + e _ 7T° 

- generates and stores unweighted events; 

- fills the histograms; 

- gives the integrated cross section within cuts. 

The program is supplemented with Gnuplot scripts for visualization of 
the histograms, produced by EKHARA. 

1.3. Related Monte Carlo programs 

There exist several MC generators for 77 physics: 



1. the code written by A. Courau [10|; used in [111, |12| for the e + e _ — > 
e + e~7r°7r°, and in [13| for the e + e~ — > e + e~7i° process; 

2. the code by F. Nguyen et al. [IH for e + e~ — > e + e~ir°n°; 



3. TREPS program written by S. Uehara [15| and used by Belle collabo 



ration 16, 17 



4. TWOGAM developed by D. M. Coffman, used by CLEO |l8| for the 
e + e~ — > e + e~P, with P = tt ,i],7]'] 

5. GGResRC used by the BaBar collaboration [19j ; 



6. GALUGA by G. A. Schuler |2fl for LEP2 physics; 



7. GaGaRes written by F. A. Berends and R. van Gulik [21(, for the study 
of resonance production in 77 interaction at LEP2 energies. 

GALUGA and GaGaRes generators describe high-energy physics and adapt- 
ing them to much lower energies is not straightforward. Other programs 
are not public. Moreover, the Equivalent Photon Approximation (EPA) is 
employed to a large extent in the majority of the generators. EPA is a use- 
ful simplification in the description of the two photon processes, when the 
accuracy requirements are not high. It leads however to some discrepancies 



with respect to the exact formulation, as shown already in [22J , especially for 
single-pseudoscalar final states. 

Summarizing, the basic requirements for a MC generator for e + e _ — > 
e + e~P studies not restricted to the region where the photons are quasi-real, 
are: 
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1. use exact formulae and exact kinematics (do not use EPA), 

2. include both s- and t-channel amplitudes and their interference, 

3. allow user-defined form factors, 

4. implement specific kinematic cuts, 

5. account for the peaking behaviour of the cross section, in order to have 
a good Monte Carlo efficiency. 

The EKHARA ver. 2, presented in this paper, fulfils all the listed above 
criteria. 

2. The generation of four-momenta for the one-pion mode 

The differential cross-section for the reaction e + e~ — > e + e~7r°, averaged 
over helicities of the initial e + e~ states, is given by 

d(x atig (e + e" -> e+e^vr ) = ~ d Lips 3 \M n o | 2 . (1) 

Here s is the initial electron-positron invariant mass squared, 1/4 is the av- 
eraging factor, 2s is the flux factor. By dLips n we denote a differential 
n— body Lorentz-invariant phase space and A4 n o stands for the matrix ele- 
ment describing the reaction e + e~ — > e + e~7r°. We generate the kinematics in 
the center-of-mass frame of the initial e + e~, with the z-axis along the initial 
positron momentum. The event is determined by five kinematic invariants 
and one azimuthal angle 0: 

da(s; t 1 ,t 2 ,s 1 ,s 2 ,(f)), 
(P1+P2) 2 , 
(Pi - <?i) 2 , 
(P2 ~ <?2) 2 , 
(Pi + 12) 2 , 

(p2 + qi) 2 - (2) 

2.1. The matrix element 

The matrix element for the reaction e + e~ — > e + e~P at the tree level 
contains the t-channel and the s-channel parts depicted in Fig. [1] (Ai n o = 
Mt + M a ). 



da(e + ( Pl ) e-(p 2 ) -> e + {q x ) e -(q 2 ) rc°(Q)) = 

s = 

ti = 

t 2 = 

51 = 

52 = 
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Q 



- ► - 




Figure 1: The f-channel (left) and the s-channel (right) diagrams for e + e 



P 



The t-channel matrix element has the following form 



M t 



F{h,t 2 



1 



-fiua/3 



fir H ti t 2 

x (v(pi) 7 M ("(92) 7^ «(P2)) 



(3) 



^0123 



1. 



The completely antisymmetric tensor ^ VOL ^ is defined by 60123 = — e 
The t-channel contribution to the cross section, is highly peaked at small t\ 
and/or £2 values. The pion two- pho ton transition form factor F(t\,t-i) is an 
important ingredient, see, e.g., [23( and it provides an additional dumping 
of the amplitude at large values of t\, t 2 . The normalization is F(0, 0) = 1, 
a ~ 1/137 is the fine structure constant and f n ~ 92.4 MeV is the pion decay 
constant. 

In Fig. |5] we demonstrate an agreement of the Monte Carlo simulation 



with the "single-tag" experimental data from CLEO [l8[ and BaBar 19 
For simulation we use the matrix element ([3]) with the form factor of the 
lowest meson dominance model with two vector resonances (LMD+V), fit- 



ted [24J to the BABAR data [19(|. A number of other expressions for the form 
factor are also available for the user of EKHARA. For example, one can use 
the formulae given by the simple vector meson dominance ansatz (rho meson 
pole), or the lowest meson dominance approach results with one vector res- 
onance multiplet (LMD) or two multiplets (LMD+V) with the generic form 



obtained from the operator product expansion considerations [23J . One may 



also use the form factor derived in a quark model, which resembles a reason- 
able agreement with the BaBar 19j data at high momentum transfer 25 . 
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Figure 2: Comparison of EKHARA Monte Carlo simulation with experimental data from 
CLEO [HI and BaBar [3] 

The matrix element for the s-channel reads 

4 io? 1 

M 8 = —r- F(s, {qi + q 2 f) — ■ ( Pl + p 2 ) a {qi + q 2 f 

In s (qi + q 2 y 

(v( Pl ) Y u(p 2 )) (u(q 2 ) Y v( qi )) . (4) 

The form factor F(s, (qi + q 2 ) 2 ) in the s-channel has to be given by the same 
analytic function as in the t-channel case. The s-channel contribution, having 
a peak at small invariant masses of the final e + e~ pair, is much smaller than 
the t-channel one when at least one of the ti, t 2 invariants is small and it 
is often missing in the existing MC codes. In EKHARA all the amplitudes 
are included (s- and t-channel as well as their interference). This guaranties 
that no relevant piece is missing even for the kinematic configurations where 
both ti and t 2 are large. 

The peaking behaviour of the s— and t— channel amplitudes is different. 
Therefore EKHARA uses, depending on the set of amplitudes included in 
the calculation (JA 8 , M.% or both), one of the three generation procedures 
described below. 

2.2. t— channel 

To generate the particle four-momenta we have adopted a very efficient 
generating procedure, which was used in the MC generator GALUGA [20]. 
Here we only sketch its layout and for further details we refer the reader to 
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201 ] . The Lorentz- invariant phase space is mapped to a unit hypercube in the 



space of the uniformly distributed random numbers G [0, 1], i = 1, . . . , 5: 

5 

dLips 3 = V* dtidt 2 dsids 2 d0 = V* J t \[ dr* , (5) 

i=l 

where Jt is the mapping Jacobian and the volume factor reads 

4vr 2 1 

v ' ■ < 6 > 

with 

/3 = 1 - Am 2 Js . (7) 

The electron mass is denoted by m e . 

The user is supposed to provide the following cuts (minima and maxima): 

• for t 2 : CUT_t2min and CUT_t2max, 

• for t x : CUT_tlmin and CUT_tlmax, 

• for positron polar angle: CUT_thlmin and CUT_thlmax, 

• for electron polar angle: CUT_th2min and CUT_th2max, 

• for positron energy: El_min and El_max, 

• for electron energy: E2_min and E2_max. 

The kinematic invariants are generated in the following sequence: 

1. calculate cuts on t% 

• t2max according to CUT_th2min and E2_min, E2_max 

• t2min according to CUT_th2max and E2_min, E2_max 

• correct t2max and t2min according to CUT_t2min and CUT_t2max 

2. generate t 2 ] t2min <t 2 < t2max 

3. calculate cuts on t\ 

• tlmax according to CUT_thlmin and El_min, El_max 

and according to generated t 2 (the lowest is taken) 
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• tlmin according to CUT_thlmax and El_min, Eljiax 

and according to generated t 2 (the highest is taken) 

• correct tlmax and tlmin according to CUT_tlmin and CUT_tlmax 

4. generate t\] tlmin < t\ < tlmax 

5. generate s\ 

6. generate s 2 

The changes of variables are the following: 

„ . / t2max\ 

t 2 = timin exp ri In , (8) 

\ t2mm J 

( tlmax \ 

t\ = tlmin exp r 2 m— — — , (9) 
\ tlmm J 

si = ^- + m 2 e +t 2 + 2m 2 e ^-, (10) 

-6 + VAsin(7r(r 4 -l/2)) 
52 = 2a ' (U) 



where 



X 1 = (u + W)(l + yi )ex P (r 3 S) , (12) 
_ s(l + /3) 2 

" ln (u + W)(l + yi )(l + y 2 ) ' {U) 

W = ^v 2 -ht 2 , (14) 

v = {ml-h-h)/2, (15) 



2/1,2 = y/l-iml/hf, (16) 

where stands for the 7r° mass. Let A 4 be the 4x4 symmetric Gram 
determinant of any four independent vectors formed out of pi, p 2 , q\, Q, 
q 2 , then its expansion in powers of s 2 determines the coefficients a, b and c: 
16 A4 = as 2 + bs 2 + c, we also use A = b 2 — 4ac. For numerically stable forms 
of A, a, b and c we refer to [201. From ti, t 2 , si, s 2 one can calculate the final 



state particles four-momenta in a frame where the x — z plane is given by the 
initial positron and final 7r° momenta. Subsequently, the 'event' is randomly 



rotated around the z— axis [20] using the r$ random number. 
For this mapping one has 

J t = 5 ln(tlmax/tlmin) In(t2max/t2min) D t = J t D t , (17) 
D t = ht 2 . (18) 
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The D t factor absorbs the peaking behaviour of the t— channel amplitude. 
The value of J t = Jt/D t for each given event is stored in the variable 
JacobianFactor. The differential cross section with the t— channel map- 
ping reads 



11 

d<V» = 4^^*11 dr * (a^I^I 2 )- (19) 



The above procedure is used when user requires to use the M.t amplitude 
only. 

2.3. s— channel 

The Lorentz-invariant phase space d Lips 3 in eq. ([T]) can be factorized in 
the following way: 

dk 2 

d Lips 3 = — — d Lips^ d Lipsg , (20) 

where k 2 has a natural meaning of the invariant mass squared of the virtual 
photon, which couples to the final e + e~ in the s— channel diagram. 
Explicit expressions for these 2-body phase spaces are 



dLi ^ = (2^v;V ' 3 + t l"'"' )2 - Fd "-- < 21 > 

1 1 fk 2 
(2tt) 2 4^P V T 



dLips£ = ———^/--^dnj, (22) 



where d is the differential solid angle in the center of mass frame of final 
e + e~, i.e. the self frame of the k 2 and dQ{ denotes the pion solid angle in 
the center of mass frame of the initial e + e~ pair. Thus, (1201) takes the form 

dLips 3 = V s dPd^dfi} , (23) 

where the volume factor reads 



V = 4-^J k - - nW<« + - *> • (24) 
2 9 vr 5 v / sFV 4 e V 4s V ' 
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The Lorentz-invariant phase space is mapped to a unit hypercube in the 
space of the uniformly distributed random numbers G [0, 1], i = 1, . . . , 5: 

5 

dLips 3 = 1/ s J s JJdr i; (25) 
i=i 

where J s is the Jacobian for a mapping from physical variables to rj. 
The following mapping is adopted: 

0, = 2vrn , (26) 
cos^ = -1 + 2r 2 ; sin 6^ = 2^/r 2 (l - r 2 ) , 
(p* f = 2nr 3 , 

cos^ = -1 + 2r 4 ; sin0* = 2^(1 - r 4 ) , 

and the invariant mass squared (k 2 ) of the final e + e~ pair is generated using 
logarithmic mapping 

e _ iml (i^iy. (27) 

From f[2"6"j) and (|2"T|) one can recover the momenta of all final state particles: 

\Q\ = ^=y/\(s,k*,ml) , (28) 

sin#j sin0. 
Q = \Q\ | sin^j cos0 
cos 0, 



and 



1 92 

-M 

92 

92* 
ft* 



i=VA(fc2,m2,mI) , (29) 

2v ft 
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The kinematic variables with asterisks are given in the center of mass 
frame of the final e + e~ pair, i.e. in the self frame of the second virtual 
photon. In order to obtain the values of the final positron and electron 
momenta in the lab frame, the q\ and q2, a Lorentz transformation of q\ and 
q\ is performed. 

For this mapping, we have 

J s = 4(2tt) 2 In ( (v/i 4 ~^ )2 ) D * = JsD s , (30) 
D s = k 2 . (31) 

The D s factor absorbs the peaking behaviour of the s— channel amplitude. 
The value of J s = J s /D s for each given event is stored in the variable 
JacobianFactor. The differential cross section with the s— channel map- 
ping reads 

11 5 

d*«* = 4^ S ^II dr * (D S J2\M\ 2 ). (32) 

i=i 

The above procedure is used when user requires to use the A4 S amplitude 
only. 

2.4- Merging s— and t— channels 

The differential cross section can be written as 

<•*«. - isn^wft)^' 2 (33) 

- \hfl dr '( AV - 3 - +BV ' J ')mfm^ 2 ^ 
= ^ J r dI ,f[ dr ^.j._^_ 5: i^ 

+ \h S?'^ Ar - v ' 3t ^^<^ m2 ' (35) 

where r, G [0,1], % = 1, . . . , 5 are the uniform random numbers, and A = 
(1— r m ), B = r m , with an a priori weight r m G (0, 1). We find empirically that 
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r m = 0.9 gives an efficient merging in order to generate events distributed 
according to \M S + M t \ 2 . The s-channel (t-channel) is used with a probability 
1 — r m (r m ). In each channel procedures described in Sections 12.21 and [273]) 
are used. 

The above procedure is followed when user requires to use the full ampli- 
tude M = M s + M t . 

2.5. Generation of the unweighted events 

Let the UB be a pre-evaluated upper bound for the Monte Carlo inte- 
grand Contrib. For generation of the unweighted events, we use the following 
accept-reject method: 

1. calculate Contrib, 

2. generate random number r accept , 

3. accept event, if Contrib > r accept UB. 

The explicit expression for Contrib depends on the mapping, which is used 
by the generator, and reads 

• s— channel 

Contrib = ^V S J S ^|-M| 2 , 

• t— channel 

Contrib = Uf l J t Y^\M\ 2 , 

• Merged s— and t— channel 

f 1 VS J* BD a +AD t T.\M\ 2 with probability (1 - r m ), 
Contrib = < 

I \ Vt Jt B S:+AD t E l-^l 2 with probability (r m ). 

The flowchart for Monte Carlo generator routine is given in Figure [31 It 
illustrates the order of phase space generation, application of cuts and the 
use of accept/reject method. In the block of the numerical stability control 
we make sure that the upper bound of the integrand is not overshooted, the 
energy and momentum conservation holds true, the final particles have on- 
shell momenta, the Gramm determinants are positively defined and that the 
generated sines and cosines of the particle angles are within [—1, 1]. 
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3. The generation of four-momenta in the two-pions mode 

The differential cross-section for the reaction e + (pi)e~ (p 2 ) — ► 
7r + (7r 1 )7r _ (7r 2 )e + (g 1 )e _ (g2), averaged over helicity states of the initial e + e _ , 
is given by 

da avg (e + e~ -)> 7i + n~e + e~) = d Lips 4 ^ \M n + n - | 2 . (36) 

The explicit expression for the matrix element Ai^^- can be found in j2| 
and its numerical evaluation is carried out in the helicity amplitudes frame- 
work. 

The multi-channel variance reduction method is used to improve effi- 
ciency of the generator and the generation is split into four channels, where 
two of them absorb peaks present in t— channel diagrams and other two take 
care of the s— channel peaks. At difference to one-pion mode only one pro- 
cedure for generation is used independently on the included contributions. 
This guaranties an efficient generation when s- and t-channel contributions 
are summed in the matrix element. Moreover, as this mode was constructed 
for generation of a background for the radiative return events, the generator 
is not optimised to run in the region, where 7*7* contributions dominate. 

In this Section we use the following definitions: 

h = qi + q 2 , Q = tti + 7r 2 , (37) 
s = (Pi + P2) 2 ,t=(q 2 - p 2 f , h = (pi - qi f. (38) 
3.1. t— channel 

For the t— channel peaks absorption, we use the following phase space 
representation 

dLips 4 (pi +p 2 ;gi,g2,vri,7r 2 ) = 

dQ' 2 dQ 2 
dLips 2 (pi + p 2 ; Q', q 2 )— — dLips 2 (Q'; Q, qi)—— dLips 2 (Q; tti, tt 2 ) (39) 

in one of the channels and an analogous one, with q\ -h- q 2 , in the other 
channel. As both channels are completely symmetric under qi -H- q 2 , we will 
describe here only changes of variables, which smoothen the distribution, 
only in one of them. For the two invariant masses (Q 2 and Q' 2 ) the following 
change of variables was performed 

Q 2 = ( v / i-2m e - v 73 ^) 2 , z = -^(^7s-2m e -2m w ) 2 (l-r Q 2) , (40) 
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Q ' 2 = v^ +m - (41) 

V ~ _ 3(g 2 + 2v^2m e ) 3 + l v 3(g 2 + 2 v ^2m e ) 3 ~ ?>{s - 2^~sm e f ' ^ 

The angles of q 2 vector are defined in the initial e + e~ center of mass (cms) 
frame with z-axis along pi and the polar angle is used to absorb the peak 
coming from the propagator of the photon exchanged in the t-channel 



3m 2 e - s + Q' 2 - 2t ^ 1 
cos 6 q2 = , t = — , 

y/l-*£\V*(s,Q»,ml) y 
1 sy/l-*£\V\s,Q",ml) 



V = + m 2 e (Q> 2 - m 2 e ) 2 ' ^ = ^ ' (42) 



where 



t = \ [Zml -s + Q' 2 -^l- ^X^(s, Q' 2 , ml) ) . (43) 

The angles of the Q vector are defined in Q' rest frame and the appropriate 
change of variables reads 

2E[Q - Q 2 - 

cos9 Q = , 0Q = 2vrr^ (44) 

2 \Pi\\Q\ 

-1 2 
2\p\\\Q\{Q 2 - 2E[Q Q - 2\p\\\Q\) + (Q 2 - 2E[Q ) 2 - Mj/^IQI 2 ^ ' 

where ^ = ^±^, Q = ^ + ^ , Ip'J = ^^11, |Q| = 
xl/2 (Q^l,rn 2 e ) and A ^ ftj c ) = a 2 + & 2 + c 2 _ 2 ( a 6 + ac + be). As we chose 

— * 

here the z-axis along the p\ (the p\ in the Q' rest frame) the vectors are 
rotated after generation to restore the general choice of the z-axis along pi 
in the e + e~ cms frame. 

Finally the angles of the positively charged pion are generated in the Q 
rest frame with flat distributions 
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cos 6 7T1 = -1 + 2r ejii , (f) wi = 



(45) 



The described change of variables transforms the phase space into a unit 
hypercube (0 < r\ < 1 , i = Q 2 , • • • , and collecting all the jacobians it 
reads 



dLips 4 (pi + p 2 ; qi, q 2 , n u vr 2 ) = P(q u q 2 )dr Q 2dr Q ,2dre q2 dr 4>q2 dr eQ dr (t>Q dre„ 1 dr^ i 

(46) 

with 



1 » i /o / ^o Ox » i /o / ^/o Ox ^ 4 771^ 



P (?i,?2) = ^ ^ , AV 2 (g' 2 ,g 2 , m 2 ) A^foQ*,^) Jl 



7T 



6(47r) 5 Q' 2 m 2 v ^ ' ^ ' e/ v ' e/ y 
g (g' 2 -m 2 ) 2 (g2_ 2g . pi)2 > /g2( v ^-2m e - 2m.) 2 

(g 2 -2^Q ) 2 -4| P 7 1 | 2 |g| 2 ' ^-2m e -y/Q* 
i i 



(g 2 + 2 A /g 2 m e ) 3 (s-2 v ^m e ) J 



(47) 



s— channel 



For the s-channel generation it is convenient to write the phase space in 
the following form 



dLips 4 G»i + p 2 ;q 1 ,q 2 ,Tr 1 ,ir 2 ) = 

dk 2 dQ 2 
d Lips 2 (pi + p 2 ; Q, h)-^ 1 d Lips 2 (fci; qi, q 2 ) -7— d Lips 2 (Q; 7Ti, 7r 2 ). (48) 

The two generation channels used here differ only in the generation of 
the electron-positron pair invariant mass k\ and the change of variables will 
be described simultaneously. The invariant mass Q 2 is generated with a flat 
distribution 

Q 2 = Ami + ((v 7 ^ - 2m e ) 2 - 4m 2 )r Q2 . (49) 
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Reflecting two leading k\ behaviors of the cross section, the two distinct 
changes of variables are done in the generation of k\: 



k\ 
Vi 



sexp(y] /3 ) , 
\n 3 (4m 2 e /s) + 



(50) 



In 3 f(l - a/QVs) 2 ) - hi 3 (4m 2 /s) 



r kl,i 



k\ = s (1 - exp(-y H )) , 
Vn : 



Q 2 (2^- v^ 2 ") 



(51) 



- ln(l - Amjs) - In I (g _ ^ I r^ n . 



The ki polar angle is used to absorb peaks coming from the electron 
propagator, while its azimuthal angle is generated with a flat distribution: 



(f) kl = 27rr 0fci , cos6» fel = 



-k\ + 2/cioPio 



2|fci||pi 



tanh 



(I) 



y = In 



kl - 2k w p 10 + 2|fci||pl| 
k\ - 2k 10 p 10 - 2|fci||pl| 



In 



kl - 2/cioPio - 2|fci||p 



kl - 2k 10 p 10 + 2 1 £4 



where fci = 2 ^ , P10 = Y> l^ 1 ' = V I ~ m e and Fil = ^75 — 

defined in the pi + P2 rest frame. 

The q\ and the 7?i angles are generated with flat distributions 



, are 



(j) qi = 2-nr^ , cos 9 Ql = -1 + 2r dqi , 0^ = 2vrr^ i , cos W1 = -1 + 2r 07ri . (53) 

After the described changes of variables are performed, the phase space 
reads (i — I or 77) 



dLips 4 (pi +p 2 ; 91,92,71-1,^2) = P s ^dr k 2dr Q 2dr eki dr ( p H dr eqi dr 4)qi dr0 ni dr^ i , (54) 
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with 

p -< = i^\/^\A r ? Al/2(s ' ' fc?)(( ^" 2m ' )2 " 4m " ) 

\h\\pi\ ( -kf + 2k 10 p 10 \ ( -kj + 2k 10 p 10 

cos 6 kl — r- — ■ h cos 9 kl 



2k loPw -k 2 i y 2|ATi||pi| ' VI 2 l*ilb~ 

-fc; + 2fciopi + 2|fc1||p 
-k 2 1 + 2k 10 p 10 -2\k 1 \\p 



, Pi , ln (zM±^o±^lMlV (55 ) 



where 



P/ = In 3 ((l - v 7 ^) 2 ) " ln 3 (4m2/ s ) , (56) 
Pu = J (izigg) ) (57) 

for the change of variables from Eq. (j50p or Eq. (1511) respectively. Again < 
Ti < 1 for i = kf, ■ ■ ■ ,(f) ni . 

3.3. Merging s— and t— channel 

The function, which approximates the peaking behavior of the matrix 
element reads 



+ — , with P s = „ ,, t , 2 . , — . (58) 



P(?l,?2) P(?2,9l) Py ' 31n^) + 



Similarly to the one pion case the cross section Eq. (156"|) is rewritten as 

1 1 



^|AW-| 2 (2 + a)Fx j 
d r 9 dr Q 2 dr Q a dr Bq2 dr^ dr dq dr^ Q dr e ^ dr^ 



da av „(e e — > n n e e 





2A 

+ / drgdr Q 2drQ / 2drg qi dr <j)qi dr eQ dr < j )Q drg^dr 4 



A 



+ / dr 9 x 



B r l 



2A UO JB 



dr 10 + / drio 



dr k 2dr Q 2drg ki dr^ dr dqi dr^ dr e ^ dr^ 



(59) 
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with A = and B = D , * . There are three channels generated with 
probabilities A, A and 1 — 2A, while the third channel branches into two 
channels generated with probabilities B and 1 — B. In the first two channels 
the change of variables described in Section 13.11 is used, while in the third 
one (branched into two channels) the changes of variables are described in 
Section 13.21 The introduced a priori weights a and b, which give the best 
efficiency of the generation at the meson factory were found to be a = 1.1 
and b = 1000. 

The flowchart for Monte Carlo generator routine is given in Figure HI 



4. Overview of the software structure 

Let us overview the directory structure of the distribution. EKHARA is 
distributed as a source code. The code of the Monte Carlo generator is lo- 
cated in the directory ekhara-routines. The main source file of EKHARA is 
ekhara.for. There are other source files in the directory ekhara-routines, 
which are automatically included: 

• the e + e~ — > e + e~7r° mode is implemented in routines_lpi . inc . f or 
and its supplementary histogramming routines are given in 
routines-histograms_lpi . inc . for; 

• the e + e~ — > e + e~7r + 7r~ mode is coded in routines_2pi . inc . f or, 
its supplementary histogramming routines are in 
routines-histograms_2pi . inc . f or and helicity-amplitude rou- 
tines are given in routines-helicity-aux . inc . f or; 

• the routines for the matrix and vector manipulations are located in 
routines-math . inc .for; 

• in routines-user . inc . for several routines are collected, which can 
be changed by a user in order to customise the operation of EKHARA; 
these include the data-card reading, reporting of events, form factor 
formulae, filling the histograms, additional phase space cuts, etc.; 

• all common blocks are included from the file common . ekhara . inc .for. 
This file contains the detailed comments on the explicit purpose of the 
most important common variables. 

The operation of the EKHARA generator requires the following steps: 
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Figure 4: Flowchart for e+e — > e+e ir + n event generation. 
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1. initialization, 

2. event generation, 

3. finalization. 

This structure allows us to build a stand-alone generator, as well as an 
interface to a separate program (e.g., a detector simulation) when EKHARA 
is called on an event-by-event basis. In the main directory of the distribu- 
tion there are examples of both uses of EKHARA. A stand-alone example is 
given in ekhara-standalone . f or. An example of an EKHARA interface on 
an event-by-event basis is given in ekhara-call-example . f or. The main 
directory of the distributed version contains a readme.txt file with a short 
description how to compile, run and test the program in the regimes, de- 
scribed above. It is suggested to use the Makefile, which is placed in the 
main directory. An example of the full set of input files and the plotting 
environment is supplied in the Env subdirectory. If one uses the distributed 
Makefile, the content of the Env subdirectory will be put into the EXE sub- 
directory together with an executable ekhara . exe (for details, see Section 
and readme.txt). In the following we assume that ekhara.exe is located 
for execution together with the input files in the EXE subdirectory. 

A source code of RANLUX random number generator^) written in C 
(ranlxd. c), together with its C-FORTRAN wrap (ranlux_f ort . c for stan- 
dard build and ranlux_f ort_vs . c for build with cl in MS Windows and xcl 
in IBM AIX) are supplied in the directory ranlux-routines. 

4-1- I/O scheme 

All the input files of EKHARA are supposed to be located in the same 
directory as the main executable, ekhara.exe. There are the following types 
of the input files: random seeds, parameter input, data-cards and histogram 
settings. An example of the full set of input files can be found in the Env 
directory. 

All the output files of EKHARA are written into . /output subdirectory. 
There are the following types of the output files: logs of execution, histograms 
and events. 



1 The unpublished double precision version of fast RANLUX 2g, |27| written by 
M. Liischer (F. James, private communication). 
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Input files 

The main input file is called input.dat. It contains all global settings: 



• number of generated events 

• mode selection (one pion or two pions in the final state) 

• histogram writing switch 

• events writing switch 

• particle masses and constants 

Another parameter, the random seed mode switch, chooses the way the seeds 
for the random number generator are handled. In the constant seed mode, 
the file seed.dat is used for every execution and it is never modified. In the 
variable seed mode, the file seed-v.dat is used and on successful completion 
of each run, a new random seed is written into seed-v.dat. The latter 
mode is convenient for a subsequent production of statistically independent 
samples. 

The channel-dependent parameters, which are supposed to be of- 
ten changed by a user are collected in "data-cards" card_lpi . dat and 
card_2pi.dat. These data-cards allow to set the total energy, types of in- 
cluded amplitudes and kinematic cuts. A detailed description can be found in 
comments within these files. In card_lpi . dat one can also use the piggFFsw 
switch in order to select the form of the two photon pion form factor. 

1 (WZWconst) constant form factor; this is not physical and should only 
be used for tests; 



2,3,4 (rho pole, LMD, LMD+V) the form factors given in 23 



5 (LMD+V new) the form factor of the lowest meson dominance model 



with two vector resonances fitted [24J to the BABAR data [19(. This 
is the recommended form factor; 



6 (quark) the form factor given in 25 



The channel-dependent histogramming settings are given in the files 
histo-settings_lpi.dat and histo-settings_2pi.dat. 
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Output: logging 

The main execution log file is the output/runf low. log. It contains the 
main information about the operation mode and status of EKHARA, this 
information is also partly written into the standard output (i.e., the console). 
At the end of a successful execution, the total cross section is reported to 
output/runf low. log and also to the standard output. 

A non-standard behaviour of the MC generator is reported into 
output /warnings . log, while the critical problems in the event generator 
operation are reported into the file output/errors . log. 

In the case of a correct operation, the output/errors . log and 
output /warnings . log should remain blank. We strongly recommend to 
keep track on this issue and report to the authors any warnings or errors. 

Output: histograms. Plotting scripts 

When histogramming is allowed through settings in the input.dat, the 
plain text files with the histogram data are saved at the end of the generator 
execution. 

• In e + e~ — > e + e~7r + 7r~ mode the file histograms_2pi . out contains 
the data for da/dQ 2 histogram. One may use the plotting script 
doplots.sh from directory histo-plotting_2pi in order to plot this 
histogram (an installed Gnuplot is required). 

• In e + e~ — > e + e~n mode there is a wide set of histograms stored in 
the files histo<Number> . <variable> .dat, where <Number> stands for 
the histogram number and <variable> is the histogramming variable 
acronym. 

One may use the plotting script do-everything. sh in the directory 
histo-plotting_lpi in order to plot all the histograms and collect 
them into a single postscript file. An installed ETjrjX system is required 
for the latter. 

One may use the plotting script doplots . sh in the directory 
tl-t2-bars_lpi in order to plot the 3D-bar graph, which shows the 
distribution in two variables: t\ and t 2 . 

As the histograms are stored as plain text files the user can use also her/his 
favourite plotting programs to visualize the histograms. 
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Output: events 

The generated four-momenta of the particles are stored in the following 
variables, which can be accessed through common blocks: 

pi initial positron, 

p2 initial electron, 

ql final positron, 

q2 final electron, 

qpion final pseudoscalar (e + e~ — > e + e~n mode), 

pil, pi2 final pseudoscalars (e + e~ — > e + e~7r + 7r~ mode). 

In the file ekhara-call-example . f or we give an example how the 
generated momenta can be used, when EKHARA works in the event-by- 
event regime. In the standalone regime we suggest to use the routine 
reportevent_lpi defined in the file routines-user . inc .for, which is called 
automatically for every accepted unweighted event (e + e~ — > e + e~7r° mode 
only). In the distributed version this routine writes the events to the file 
output/events . out when WriteEvents flag is on. One can modify this rou- 
tine in order to accommodate the way how the generated events are collected. 

4-2. Selected procedures 

The top-level interface to the Monte Carlo generator is provided by the 
routine 

EKHARA (i) i = -1: initialize, 

i = 0: generate event (s), 
i = 1: finalize. 



Only this routine should be called from an external program, 
when you use EKHARA in the event-by-event regime, see example 
ekhara-call-example .for. 

In order to describe briefly the "internal" structure of EKHARA, we list 
several important routines. 

EKHARA_INIT_read reading the input files and datacards, 

EKHARA_INIT_set initialization of the MC loop and mappings, 

EKHARA_RUN MC loop execution, 

EKHARA_FIN MC finalization, saving the results. 



25 



e + e — > e + e ir° procedures 



mc_loop_lpi 
EvalUpperBound 

phasespace_lpi 
event select ion_lpi 
m_el_lpi 
tellSIGMA_lpi 



Monte Carlo loop (see the flowchart in Figure [3]) , 
evaluation of the upper bound for the Monte Carlo 
integrand, 

a wrap for the phase space generation routines, 
kinematic cuts, 

calculation of the matrix element for e + e~ — > e + e~7i°, 
reports the total cross section. 



e + e — > e + e 7r + 7r procedures 



mc_loop_2pi 

event select ion_2pi 

phspl 

phsp2 

phsp3 

matrixelmt 



Monte Carlo loop (see the flowchart in Figure H|) , 
kinematic cuts, 

phase space generation routine (branch 1), 
phase space generation routine (branch 2), 
phase space generation routine (branch 3), 
calculation of the matrix element for e + e~ 
e + e~7t + TT~ , 

reports the total cross section. 



tellSIGMA_2pi 

For details see the comments in the source code. 



5. Compilation instructions 

Being distributed as a source code the program does not require installa- 
tion, but a compilation and linking is needed. EKHARA does not need any 
specific external libraries. In order to compile the program, a user may run 
any OS (UNIX, Linux, MS Windows, etc) with correctly installed 

• FORTRAN 77 compiler with support of quadruple precision, 

• C compiler. 

The program was tested on the following platforms: 

• Linux (Ubuntu 8.04.3) 
GNU C Compiler (gcc) 4.2.4 

Intel(R) FORTRAN Compiler (ifort) 11.0.20080930 
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Description 



Linux Windows IBM AIX 



Default: Standalone MC gen- default default-vs def ault-ibm 
erator, sample input files and 
histogramming routines. Ev- 
erything put into EXE directory 

Compile everything: default, all all-vs all-ibm 

ranlux-testing program and 
seed-production programs 

Testrun: Compile everything test test-vs test-ibm 

and execute the testrun scripts 

Remove the redundant and clean clean-vs clean-ibm 

temporary files 



Table 7: The main Makefile targets 



• MS Windows (XP SP3) 

MS Visual Studio 2008 (nmake, cl) 
Intel(R) FORTRAN Compiler (ifort) 11. 

The program distribution contains the Makefile, with targets for Linux, 
IBM AIX and Windows environments. The Makefile is annotated, in order 
to help a user to tune it up for the own requirements. The main Makefile 
targets are listed in Table [3 

For example, in Linux, a simple way to compile a program is to issue 
make default, being in the directory where the Makefile is located. This 
will produce ekhara . exe (main program executable) and copy it into the sub- 
directory EXE, together with the contents of Env sub-directory. The latter 
contains the set of sample input files and histogram plotting scripts. We 
provide a full set of necessary input files in the distribution package. It 
is advised to execute ekhara.exe in the directory EXE, where it is placed 
by default. Every time one does make default, the input files in EXE are 
replaced with the sample ones from Env. 
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In order to produce only an object file with the EKHARA MC generator, 
one can use, for example 

ifort -c ekhara-routines/ekhara. f or -o ./ekhara.o 

EKHARA needs a random seed for operation. Different random seeds can 
be obtained by using a Makefile target seecLprod-if ort. It produces an 
exacutable program seed_prod.exe, which generates a set of random seeds. 

6. Test run description 

It is recommended to test the random number generator on a given 
machine, before using EKHARA. It is also important to check whether 
EKHARA can function properly on a given operational system and that 
there are no critical bugs due to the compiler. We provide a testrun package 
for these purposes. 

It is suggested to use the Makefile targets test, test-vs or test-ibm 
depending on your environment (Linux, Windows and IBM AIX, correspond- 
ingly). This will automatically prepare and execute the following two test 
steps. 

The first step of the testrun is the random number generator control. The 
source file testlxf .for contains the ranlux test routines. 

The second step is the verification if the user-compiled EKHARA can re- 
produce the set of results, created by a well-tested copy of EKHARA in vari- 
ous modes. The testrun environment contains directories test and test-vs 
with precalculated data for a comparison, along with the random seed and 
input files for each mode. The script test . sh is responsible for the exe- 
cution of a user-compiled ekhara.exe in all the control modes and for the 
comparison of the output. 

Please read carefully the output of the testrun execution in your console 
and be sure there are no warnings and/or error messages. 

7. Customization of the source code by a user 

We leave for a user an option to customise the generator to her/his needs 
by editing the source code file ekhara-routines/routines-user . inc . f or. 
Notice, we always use explicit declaration of identifiers and the implicit 
none statement is written down in each routine. 

In the file ekhara-routines/routines-user . inc . for one can change 
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• the data-card reading (routines reacLcarcLlpi and read_card_2pi), 

• the form-factor formula (routine piggFF), 

• the events reporting (routine reportevent_lpi), 

• histogramming (routines histo_event_lpi and histo_event_2pi), 

• additional kinematic cuts (routines ExtraCuts_lpi and 
ExtraCuts_2pi). 

8. Validation of the generator 

In Fig. [2] we demonstrated an agreement of the Monte Carlo simulation 
of e + e~ — > e + e~7r° in the "single-tag" mode with the experimental data from 
CLEO 18 1 and BaBar [l9[. We conclude that the matrix element is well 



under control and the applied pion transition form factor (LMD+V) [24| is 
in agreement with data. The procedure of the phase space generation was 
validated by means of the high statistics phase space volume calculation in 
EKHARA and comparison of the result with that from the independent dedi- 
cated numerical calculation. The volume was also compared to that obtained 
by GALUGA generator [2(|. We have also verified that our numerical results 
for the three-body phase space volume in the limit of massless ir° reproduce 
well those of the analytic expression. 

The e + e~ — > e + e~7r + 7r~ mode is validated by means of the reproduction 
of the results from the previous version of EKHARA. In [H, 0, 0] the tests of 
this part are presented in detail. 

9. Summary 

An update (version 2.0) of the Monte Carlo event generator EKHARA is 
presented. It generates processes e + e _ — > e + e~n° and e + e~ — > e + e~7r + 7r~ . 
The newly added channel (e + e~ — > e + e~7r°) is important for 7*7* physics and 
can be used for the pion transition form factor studies at meson factories. 
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